Thermal Sensitivity Perspective
Pieter Arnold
2024-12-04
library(tidyverse)
library(ggpubr)
library(raster)
library(ncdf4)
library(ozmaps)
###
# Functions needed for Thermal Sensitivity Perspective Code
###
###
# Find the intercept from GLM
findInt <- function(model, value) {
function(x) {
predict(model, data.frame(log10timeh = x), type = "response") - value
}
}
###
# Calculate the accumulated damage based on time in minutes
# Function Source: Faber, Ørsted, Ehlers (2024) https://doi.org/10.1093/jxb/erae096
Acc_function = function(x, y, intercept, slope) {
output = data.frame(z = NA, cumsum = NA)
last_row = nrow(data.frame(x = x, y = y))
x = c(x, x[last_row])
y = c(y, y[last_row])
for (i in 1:length(x) - 1) {
#print(paste0("round nr. ", i))
acc = max(0, min(100, (100 * (x[i + 1] - x[i])) / (10^(slope * max(y[i + 1], y[i]) + intercept))))
output[i, 1] = acc
}
output$cumsum = cumsum(output$z)
z = output$cumsum
z = append(z, 0)
z[z > 100] = 100 # Set values exceeding 100% to 100%
z = z[-length(z)]
return(z)
}
###
# Repair function ArrFunc5 by Michael Kearney
ArrFunc5 <- function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref) {
exp(TA * (1/TREF-1/(273.15+x)))/(1+exp(TAL * (1/(273.15+x)-1/TL)) + exp(TAH * (1/TH-1/(273.15+x)))) * kdot_ref
}
###
# Function to estimate thermal tolerance landscape from static assays
# General procedure
# Step 1: Calculate CTmax and z from TDT curve
# Step 2: Calculate average log10 time and Ta (mean x and y for interpolation purposes)
# Step 3: Interpolating survival probabilities to make them comparable across treatments
# Step 4: Overlap all survival curves into a single one by shifting each curve to mean x and y employing z
# Step 5: Build expected survival curve with mean x and y pooling all data
# Step 6: Expand expected curve to multiple Ta (with 0.1ºC difference for predictive purposes)
tolerance.landscape <- function(ta,time){
data <- data.frame(ta,time)
data <- data[order(data$ta,data$time),]
# Step 1: Calculate CTmax and z from TDT curve
ta <- as.numeric(levels(as.factor(data$ta)))
model <- lm(log10(data$time) ~ data$ta); summary(model)
ctmax <- -coef(model)[1]/coef(model)[2]
z <- -1/coef(model)[2]
# Step 2: Calculate average log10 time and Ta (mean x and y for interpolation purposes)
time.mn <- mean(log10(data$time))
ta.mn <- mean(data$ta)
# Step 3: Interpolating survival probabilities to make them comparable across treatments
time.interpol <- matrix(NA,1001,length(ta))
for(i in 1:length(ta)){
time <- c(0,sort(data$time[data$ta==ta[i]]))
p <- seq(0,100,length.out = length(time))
interp <- approx(p,time,n = 1001)
time.interpol[,i] <- interp$y
}
# Step 4: Overlap all survival curves into a single one by shifting each curve to mean x and y employing z
# Step 5: Build expected survival curve with median survival time for each survival probability
shift <- (10^((ta - ta.mn)/z))
time.interpol.shift <- t(t(time.interpol)*shift)[-1,]
surv.pred <- 10^apply(log10(time.interpol.shift),1,median)
# Step 6: Expand predicted survival curves to measured Ta (matrix m arranged from lower to higher ta)
# Step 7: Obtain predicted values comparable to each empirical measurement
m <- surv.pred*matrix ((10^((ta.mn - rep(ta, each = 1000))/z)), nrow = 1000)
out <-0
for(i in 1:length(ta)){
time <- c(0,data$time[data$ta==ta[i]])
p <- seq(0,100,length.out = length(time))
out <- c(out,approx(seq(0,100,length.out = 1000),m[,i],xout=p[-1])$y)
}
data$time.pred <- out[-1]
colnames(m) <- paste("time.at",ta,sep=".")
m <- cbind(surv.prob=seq(1,0.001,-0.001),m)
par(mfrow=c(1,2),mar=c(4.5,4,1,1),cex.axis=1.1)
plot(-10,-10,las=1,xlab="Time (min)",ylab="Survival (%)",col="white",xaxs="i",yaxs="i",xlim=c(0,max(data$time)*1.05),ylim=c(0,105))
for(i in 1:length(ta)){
time <- c(0,sort(data$time[data$ta==ta[i]])); p <- seq(100,0,length.out = length(time))
points(time,p,pch=21,bg="black",cex=0.5)
time <- c(0,sort(data$time.pred[data$ta==ta[i]]))
points(m[,i+1],100*m[,1],type="l",lty=2)}
segments(max(data$time)*0.7,90,max(data$time)*0.8,90,lty=2)
text(max(data$time)*0.82,90,"fitted",adj=c(0,0.5))
plot(log10(data$time.pred),log10(data$time),pch=21,bg="black",cex=0.5,lwd=0.7,las=1,xlab="Fitted Log10 time",ylab="Measured Log10 time")
abline(0,1,lty=2)
rsq <- round(summary(lm(log10(data$time) ~ log10(data$time.pred)))$r.square,3)
text(min(log10(data$time.pred)),max(log10(data$time)),substitute("r"^2*" = "*rsq),adj=c(0,1))
list(ctmax = as.numeric(ctmax), z = as.numeric(z), ta.mn = ta.mn, S = data.frame(surv=seq(0.999,0,-0.001),time=surv.pred),
time.obs.pred=cbind(data$time,data$time.pred), rsq = rsq)
}
###
# Function dynamic.landscape1 by Michael Kearney and Enrico Rezende
dynamic.landscape1 <- function(ta,tolerance.landscape,kdot_ref){
ArrFunc5 <- function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref){
exp(TA*(1/TREF-1/(273.15+x)))/(1+exp(TAL*(1/(273.15+x)-1/TL))+exp(TAH*(1/TH-1/(273.15+x))))*kdot_ref
}
TA <- 1000
TL <- 8.5 + 273.15
TH <- 25.5 + 273.15
TAL <- 100000
TAH <- 100000
TREF <- 20 + 273.15
Tbs <- seq(0, 50)
TempCorr <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, 0.001)
surv <- tolerance.landscape$S[,2] + seq(0,0.0000001,length.out=1000)
ta.mn <- tolerance.landscape$ta.mn
z <- tolerance.landscape$z
shift <- 10^((ta.mn - ta)/z)
time.rel <- 0
alive <- 100
survival <- matrix(NA, length(ta))
repair <- ArrFunc5(ta, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
for(i in 1:length(ta)){
if(!is.na(alive)) {
alive <- try(approx(c(0,shift[i]*surv),seq(100,0,length.out = length(c(0,surv))),
xout = time.rel[i] + 1)$y,silent=TRUE)
alive <- alive + ArrFunc5(ta[i], TA, TAL, TAH, TL, TH, TREF, kdot_ref)
alive <- ifelse(alive > 100,100,alive)
survival[i] <- alive
time.rel <- try(c(time.rel,approx(seq(100,0,length.out = length(c(0,surv))),c(0,shift[i + 1]*surv),
xout = alive)$y),silent=TRUE)}
else{
alive <- 100
survival[i] <- 0}}
out <- data.frame(cbind(ta=ta[1:(length(survival)-1)],time=(1:length(ta))[1:(length(survival)-1)],
alive=survival[1:(length(survival)-1)]))
par(mar=c(4,4,1,1),mfrow=c(1,2))
plot(1:length(ta),ta,type="l",xlim=c(0,length(ta)),ylim=c(min(ta),max(ta)),col="black",lwd=1.5,las=1,
xlab = "Time (min)", ylab = "Temperature (ºC)")
plot(out$time,out$alive,type="l",xlim=c(0,length(ta)),ylim=c(0,100),col="black",lwd=1.5,las=1,
xlab = "Time (min)", ylab = "Survival (%)")
list(time = out$time,ta = out$ta, survival = out$alive)}
# Add in a dynamic repair function within by Enrico Rezende, Michael Kearney, Pieter Arnold
dynamic.landscape2 <-
function(ta,tolerance.landscape,kdot_ref){
ArrFunc5 <- function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref){
exp(TA*(1/TREF-1/(273.15+x)))/(1+exp(TAL*(1/(273.15+x)-1/TL))+exp(TAH*(1/TH-1/(273.15+x))))*kdot_ref
}
TA <- 1000
TL <- 8.5 + 273.15
TH <- 25.5 + 273.15
TAL <- 100000
TAH <- 100000
TREF <- 20 + 273.15
Tbs <- seq(0, 50)
TempCorr <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, 0.001)
surv <- tolerance.landscape$S[,2] + seq(0,0.0000001,length.out=1000)
ta.mn <- tolerance.landscape$ta.mn
z <- tolerance.landscape$z
shift <- 10^((ta.mn - ta)/z)
time.rel <- 0
alive <- 100
survival <- matrix(NA, length(ta))
kdot <- matrix(NA, length(ta))
repair <- ArrFunc5(ta, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
for(i in 1:length(ta)){
if(!is.na(alive)) {
alive <- try(approx(c(0,shift[i]*surv),seq(100,0,length.out = length(c(0,surv))),
xout = time.rel[i] + 1)$y,silent=TRUE)
# modify kdot_ref based on alive. If alive > 99, kdot_ref stays the same. If alive < 99, kdot_ref reduces as a function of alive
kdot_ref1 <- ifelse(alive < 99, kdot_ref*(kdot_ref^(3/(alive))), kdot_ref)
alive <- alive + ArrFunc5(ta[i], TA, TAL, TAH, TL, TH, TREF, kdot_ref1)
alive <- ifelse(alive > 100,100,alive)
survival[i] <- alive
kdot[i] <- kdot_ref1
time.rel <- try(c(time.rel,approx(seq(100,0,length.out = length(c(0,surv))),c(0,shift[i + 1]*surv),
xout = alive)$y),silent=TRUE)}
else{
alive <- 100
survival[i] <- 0
kdot[i] <- kdot_ref }}
out <- data.frame(cbind(ta=ta[1:(length(survival)-1)],time=(1:length(ta))[1:(length(survival)-1)],
alive=survival[1:(length(survival)-1)],kdot=kdot[1:(length(survival)-1)]))
par(mar=c(4,4,1,1),mfrow=c(1,2))
plot(1:length(ta),ta,type="l",xlim=c(0,length(ta)),ylim=c(min(ta),max(ta)),col="black",lwd=1.5,las=1,
xlab = "Time (min)", ylab = "Temperature (ºC)")
plot(out$time,out$alive,type="l",xlim=c(0,length(ta)),ylim=c(0,100),col="black",lwd=1.5,las=1,
xlab = "Time (min)", ylab = "Survival (%)")
list(time = out$time,ta = out$ta, survival = out$alive, kdot = out$kdot)}
###
# Function to extract data from netcdf files by Michael Kearney
extract.nc <- function(file, loc){
nc <- nc_open(file)
lon <- ncvar_get(nc, "longitude") # extract longitude values
lat <- ncvar_get(nc, "latitude") # extract latitude values
# find closest pixel to site requested
flat <- match(abs(lat-loc[2])<.3,1)
latindex <- which(flat %in% 1)
flon <- match(abs(lon-loc[1])<.3,1)
lonindex <- which(flon %in% 1)
start <- c(1,latindex,lonindex) # start location in file for extraction
count <- c(-1, 1, 1) # extract all data (-1)
data <- as.numeric(ncvar_get(nc, varid = attributes(nc$var)$names, start = start, count)) # get the data
# extract dates
dates <- as.POSIXct(ncvar_get(nc, "time"), tz = "Etc/GMT+10", origin = "1970-01-01")
nc_close(nc)
result <- cbind(dates, as.data.frame(data)) # final output
}
#source("Thermal sensitivity functions.R")
cbp2 <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")weeddata <- read.csv("weedseeddata_all.csv")
# Data extracted from Dahlquist et al. (2007) Figures using metaDigitise
# https://doi.org/10.1614/WS-04-178.1
weeddata$mortality_prop <- weeddata$mortality_percent/100
weeddata$log10timemin <- log10(weeddata$time_min)
weeddata$log10timeh <- log10(weeddata$time_h)
weeddata_sp1 <- subset(weeddata, species=="London rocket")
# Load in extracted dataset that contains LT80 values
weedLT80 <- read.csv("Weed_LT80.csv")
weedLT80$log10timemin <- log10(weedLT80$time_min)
weedLT80$log10timeh <- log10(weedLT80$time_h)model50 <- glm(mortality_prop ~ log10timeh,
data = weeddata_sp1[weeddata_sp1$temp_treatment=="50",], family = quasibinomial)
model46 <- glm(mortality_prop ~ log10timeh,
data = weeddata_sp1[weeddata_sp1$temp_treatment=="46",], family = quasibinomial)
model42 <- glm(mortality_prop ~ log10timeh,
data = weeddata_sp1[weeddata_sp1$temp_treatment=="42",], family = quasibinomial)# Fit overall and species by species TDT models
overallLT80_mod <- lm(temp_treatment ~ log10timeh, data = weedLT80)
LT80mod <- lm(temp_treatment ~ log10timeh, data = weedLT80)
LT80mod_sp3 <- lm(temp_treatment ~ log10timeh, data = subset(weedLT80, weedLT80$Species2=="Black nightshade"))
LT80mod_sp5 <- lm(temp_treatment ~ log10timeh, data = subset(weedLT80, weedLT80$Species2=="London rocket"))
LT80mod_sp6 <- lm(temp_treatment ~ log10timeh, data = subset(weedLT80, weedLT80$Species2=="Tumble pigweed"))
# Fit the alternative model form: time ~ temp to apply accumulated damage model
LT80mod_Orsted <- lm(log10timeh ~ temp_treatment, data = weedLT80)
LT80mod_Orsted1 <- lm(log10timeh ~ temp_treatment, data = subset(weedLT80, weedLT80$Species2=="Black nightshade"))
LT80mod_Orsted2 <- lm(log10timeh ~ temp_treatment, data = subset(weedLT80, weedLT80$Species2=="London rocket"))
LT80mod_Orsted3 <- lm(log10timeh ~ temp_treatment, data = subset(weedLT80, weedLT80$Species2=="Tumble pigweed"))
intercept <- as.numeric(LT80mod_Orsted1$coefficients[1])
slope <- as.numeric(LT80mod_Orsted1$coefficients[2])
colfunc <- colorRampPalette(c("palegreen","tan4"))
facet.labs <- c("Solanum nigrum", "Sisymbrium irio", "Amaranthus albus")
names(facet.labs) <- c("Acc_damage_sp1", "Acc_damage_sp2", "Acc_damage_sp3")
# Simulate 50ËšC temperature ramp and hold treatment with thermal variability
#st1 <- c(seq(30,39,0.5),rep(40,20),seq(41,45,0.5),rep(46,20),seq(47,51,0.5),rep(52,100))
st1 <- c(seq(20,50,2.5),rep(50,100))
specific_temperatures1 <- st1 + rnorm(st1,0,1)
Temperature <- specific_temperatures1
Time <- specific_temperatures1
data <- data.frame(Time,Temperature)
data <- data %>% filter(row_number() %% 2 != 1)
data$Time <- seq(1, nrow(data))
#data$Time <- as.numeric(data$Time)/2
data$Temperature <- as.numeric(data$Temperature)
data <- data[1:100,]
# Calculate accumulated damage
data$Acc_damage_mean <- Acc_function(data$Time,data$Temperature,
as.numeric(LT80mod_Orsted$coefficients[1]), # intercept
as.numeric(LT80mod_Orsted$coefficients[2])) # slope
data$Acc_damage_mean <- as.numeric(data$Acc_damage_mean)
data$Acc_damage_sp1 <- Acc_function(data$Time,data$Temperature,
as.numeric(LT80mod_Orsted1$coefficients[1]), # intercept
as.numeric(LT80mod_Orsted1$coefficients[2])) # slope
data$Acc_damage_sp2 <- Acc_function(data$Time,data$Temperature,
as.numeric(LT80mod_Orsted2$coefficients[1]), # intercept
as.numeric(LT80mod_Orsted2$coefficients[2])) # slope
data$Acc_damage_sp3 <- Acc_function(data$Time,data$Temperature,
as.numeric(LT80mod_Orsted3$coefficients[1]), # intercept
as.numeric(LT80mod_Orsted3$coefficients[2])) # slope
datalong <- data %>% pivot_longer(cols = c(Acc_damage_sp1, Acc_damage_sp2, Acc_damage_sp3),
names_to = "Acc_damage")
datalong <- na.omit(datalong)preddata <- matrix(data = NA, 51, 3)
preddata[,2] <- seq(0,5,0.1)
preddata <- as.data.frame(preddata)
colnames(preddata) <- c("test", "log10timeh","out")
# CTmax1h temperature (at log10(0))
predict(LT80mod, as.data.frame(preddata))[1]## 1
## 57.13306LT80mod$coefficients[1] # CTmax1h equivalent (just checking the predict works)## (Intercept)
## 57.13306LT80mod_sp3$coefficients[1]## (Intercept)
## 60.09522LT80mod_sp5$coefficients[1]## (Intercept)
## 53.47222LT80mod_sp6$coefficients[1]## (Intercept)
## 59.32692# Find time for LT80 in log10 hours
uniroot(findInt(model50, 0.8), range(weeddata_sp1$log10timeh))$root## [1] 0.60383uniroot(findInt(model46, 0.8), range(weeddata_sp1$log10timeh))$root## [1] 1.305468uniroot(findInt(model42, 0.8), range(weeddata_sp1$log10timeh))$root## [1] 1.929779# Find time for LT80 in hours
10^(uniroot(findInt(model50, 0.8), range(weeddata_sp1$log10timeh))$root)## [1] 4.01633610^(uniroot(findInt(model46, 0.8), range(weeddata_sp1$log10timeh))$root)## [1] 20.205410^(uniroot(findInt(model42, 0.8), range(weeddata_sp1$log10timeh))$root)## [1] 85.07045# Hours at treatment of 50ËšC to reach LT80
10^(uniroot(findInt(LT80mod, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # Average## [1] 19.6840910^(uniroot(findInt(LT80mod_sp3, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # Black nightshade## [1] 35.0792910^(uniroot(findInt(LT80mod_sp5, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # London rocket## [1] 3.71918810^(uniroot(findInt(LT80mod_sp6, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # Tumble pigweed## [1] 61.95814# LT80 for london rocket seeds
fig1A <-
ggplot(data = weeddata_sp1, aes(y = mortality_prop, x = log10timeh)) +
geom_point(aes(colour = as.factor(temp_treatment)), alpha = 0.4) +
geom_smooth(aes(colour = as.factor(temp_treatment), fill = as.factor(temp_treatment)),
method = "glm", formula = y ~ x,
method.args = list(family = "quasibinomial"), se = T, fullrange = T) +
geom_hline(yintercept = 0.8, lty = 2) +
geom_segment(x = 1.929779, xend = 1.929779, y = 0.8, yend = -Inf, lty = 2) + # 42
geom_segment(x = 1.305468, xend = 1.305468, y = 0.8, yend = -Inf, lty = 2) + # 46
geom_segment(x = 0.60383, xend = 0.60383, y = 0.8, yend = -Inf, lty = 2) + # 50
scale_colour_manual(values = c("goldenrod1", "darkorange2", "red3"), "Temperature (°C)") +
scale_fill_manual(values = c("goldenrod1", "darkorange2", "red3"), "Temperature (°C)") +
annotate("text", x = 0.4, y = 1.1, label = bquote(italic(Sisymbrium~irio))) +
# annotate("text", x = 0.4, y = 0.58, label = "50°C") +
# annotate("text", x = 1.1, y = 0.58, label = "46°C") +
# annotate("text", x = 1.7, y = 0.58, label = "42°C") +
annotate("text", x = 0, y = 0.75, label = bquote(italic(LT)[80])) +
annotate("text", x = 2.5, y = 0.75, label = bquote(italic(LT)[80])) +
labs(y = "Seed mortality (proportion)", x = bquote(log[10]~Time~(h))) +
scale_x_continuous(limits = c(-0.4,3)) +
scale_y_continuous(limits = c(0,1.2), breaks = c(0,0.2,0.4,0.6,0.8,1)) +
theme_classic() + theme(legend.position = "bottom", legend.box = "vertical")
#fig1A
# TDT models fit to 3 weed species and overall, CTmax' and z
fig1B <-
ggplot(data = weedLT80, aes(y = temp_treatment, x = log10(time_h),
group = Species, fill = Species, shape = Species)) +
geom_abline(intercept = LT80mod_sp3$coefficients[1], slope = LT80mod_sp3$coefficients[2],
colour = cbp2[6], lty = 1, linewidth = 1) + # black_nightshade, solanum nigrum
geom_abline(intercept = LT80mod_sp5$coefficients[1], slope = LT80mod_sp5$coefficients[2],
colour = cbp2[4], lty = 1, linewidth = 1) + # london_rocket, sisymbrium irio
geom_abline(intercept = LT80mod_sp6$coefficients[1], slope = LT80mod_sp6$coefficients[2],
colour = cbp2[8], lty = 1, linewidth = 1) + # tumble_pigweed, amaranthus albus
#geom_abline(intercept = overallLT80_mod$coefficients[1], slope = overallLT80_mod$coefficients[2],
# colour = cbp2[1], lty = 1, linewidth = 1) + # overall
geom_point(size = 3) +
geom_vline(xintercept = 0, lty = 2) +
scale_fill_manual(values = c(cbp2[6], cbp2[4], cbp2[8], cbp2[1]), name = "",
labels = c(bquote(italic(Solanum~nigrum)),
bquote(italic(Sisymbrium~irio)),
bquote(italic(Amaranthus~albus)))) +
scale_shape_manual(values = c(22,23,24), name = "",
labels = c(bquote(italic(Solanum~nigrum)),
bquote(italic(Sisymbrium~irio)),
bquote(italic(Amaranthus~albus)))) +
annotate("text", x = 0.2, y = 60, label = bquote(italic(CT[max1h])~"="~intercept), hjust = 0) +
annotate("text", x = 1.6, y = 54, label = bquote(italic(z)~"="~slope)) +
annotate("text", x = 0.05, y = 63, label = bquote(italic(LT)[80]~TDT~curve), hjust = 0) +
annotate("text", x = 0.05, y = 42.5, label = bquote(italic(CT[max1h])~"="~"60.1°C"), colour = cbp2[6], hjust = 0) +
annotate("text", x = 0.05, y = 41, label = bquote(italic(CT[max1h])~"="~"53.5°C"), colour = cbp2[4], hjust = 0) +
annotate("text", x = 0.05, y = 39.5, label = bquote(italic(CT[max1h])~"="~"59.3°C"), colour = cbp2[8], hjust = 0) +
annotate("text", x = 0.05, y = 38, label = bquote(italic(z)==-6.534), colour = cbp2[6], hjust = 0) +
annotate("text", x = 0.05, y = 36.5, label = bquote(italic(z)==-6.087), colour = cbp2[4], hjust = 0) +
annotate("text", x = 0.05, y = 35, label = bquote(italic(z)==-5.205), colour = cbp2[8], hjust = 0) +
scale_x_continuous(limits = c(-0.2,3)) +
scale_y_continuous(limits = c(35,65), breaks = c(35,40,45,50,55,60,65)) +
labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
theme_classic() + theme(legend.position = "bottom")
#fig1B
# Damage accumulation model applied to 3 species over temperature-time treatment simulation
fig1C <-
ggplot(datalong, aes(x = Time, y = Temperature)) +
facet_wrap(~Acc_damage, nrow = 3, labeller = labeller(Acc_damage = facet.labs)) +
geom_line() +
geom_point(data = datalong, aes(x = Time, y = Temperature, fill = value),
shape = 21, size = 1.5, alpha = 1) +
scale_fill_gradientn("Damage (%)", colours = colfunc(21), limits = c(0,100)) +
labs(x = "Time (h)", y = bquote(Temperature~(degree*C))) +
scale_x_continuous(limits = c(0,60), breaks = c(0,10,20,30,40,50)) +
ylim(20,55) +
theme_classic() +
theme(legend.position = "right", strip.background = element_blank(),
strip.text = element_text(face = "italic"))
#fig1C
# Damage accumulation model applied to 3 species over temperature-time treatment simulation
fig1D <-
ggplot(data, aes(x = Time, y = Acc_damage_mean)) +
geom_line(data = data, aes(x = Time, y = Acc_damage_sp1, colour = cbp2[6]), linewidth = 1) +
geom_line(data = data, aes(x = Time, y = Acc_damage_sp2, colour = cbp2[4]), linewidth = 1) +
geom_line(data = data, aes(x = Time, y = Acc_damage_sp3, colour = cbp2[8]), linewidth = 1) +
#geom_point(data = data, aes(x = Time, y = Acc_damage_sp1, fill = cbp2[6]), shape = 22, size = 2, alpha = 0.5) +
#geom_point(data = data, aes(x = Time, y = Acc_damage_sp2, fill = cbp2[4]), shape = 23, size = 2, alpha = 0.5) +
#geom_point(data = data, aes(x = Time, y = Acc_damage_sp3, fill = cbp2[8]), shape = 24, size = 2, alpha = 0.5) +
scale_fill_manual(values = c(cbp2[6], cbp2[4], cbp2[8]), name = "",
labels = c(bquote(italic(Solanum~nigrum)),
bquote(italic(Sisymbrium~irio)),
bquote(italic(Amaranthus~albus)))) +
scale_colour_manual(values = c(cbp2[6], cbp2[4], cbp2[8]), name = "",
labels = c(bquote(italic(Solanum~nigrum)),
bquote(italic(Sisymbrium~irio)),
bquote(italic(Amaranthus~albus)))) +
scale_y_continuous(limits = c(0,110), breaks = c(0,25,50,75,100)) +
scale_x_continuous(limits = c(0,60), breaks = c(0,10,20,30,40,50)) +
geom_hline(yintercept = 100, linetype = "dashed") +
labs(x = "Time (h)", y = expression(paste("Accumulated damage"," (%)"))) +
theme_classic() + theme(legend.position = "none")
#fig1D
ggarrange(fig1A, fig1B, fig1C, fig1D,
nrow = 2, ncol = 2, align = "v", labels = c("(a)", "(b)", "(c)", "(d)"))
#ggsave("Fig1_R1.pdf", height = 8.5, width = 10)set.seed(21)
# Set Arrhenius model and repair rate parameters
TA <- 14065 # Arrhenius temperature, K
TL <- 10.5 + 273.15 # Arrhenius temperature lower threshold, K
TH <- 28.5 + 273.15 # Arrhenius temperature upper threshold, K
TAL <- 50000 # Arrhenius temperature lower, K
TAH <- 100000 # Arrhenius temperature upper, K
TREF <- 20 + 273.15 # reference temperature, K
Tbs <- seq(0, 50) # sequence of temperatures over which to model, deg C
kdot_ref <- 0.02 # repair rate % per minute at ref temp 20 deg C
# Load ectotherm model data and subset to example four week period
par(mfrow = c(1,1))
ta <- read.csv("ectotherm.csv")[3006 + 1:(24*28),] # subset ectotherm - derived from NicheMapR
ta <- ta$TC + 0.3 # Body temperature
plot(ta, type = "l")
time.min <- spline(ta, n = length(ta-1)*60)$x
ta.min <- spline(ta,n = length(ta-1)*60)$y
ta <- ta.min
ind <- 1000
ctmax <- 43.8 # CTmax1h
z <- -2.4 # z (slope)
static <- rep(c(round(ctmax,2),round(ctmax+z,2), round(ctmax+z*2,2)), each = ind)
time <- abs(c(rnorm(ind,1,1/4), rnorm(ind,10,10/4), rnorm(ind,100,100/4)))
tl <- tolerance.landscape(static, time)
# Fit modified dynamic.landscape function that includes the Arrhenius repair function
# Set repair coefficients
kdot0 <- 0 # none
kdot1 <- 0.007 # low
kdot2 <- 0.0111 # moderate
kdot3 <- 0.02 # high
# No repair
dl0 <- dynamic.landscape2(ta,tl,kdot0)
# Three levels of repair rates
dl1 <- dynamic.landscape2(ta,tl,kdot1)
dl2 <- dynamic.landscape2(ta,tl,kdot2)
dl3 <- dynamic.landscape2(ta,tl,kdot3)
# Temperature-time model sequence
time2 <- log10(seq(1,330,1))
ht <- 43.8 + (-2.4*time2)
sim <- data.frame(cbind(time2,ht))
Tbs <- seq(5, 45, 0.2)
# Generate Arrhenius functions for repair rates
TempCorr0 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot0)
TempCorr1 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot1)
TempCorr2 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot2)
TempCorr3 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot3)
# Combine dataset
damage <- c(seq(0.000001,0.005,0.00005),(seq(0.1,3,0.01)^2.3))/15
damage <- damage[1:201]
range(damage)## [1] 6.666667e-08 8.300646e-02TempCorrDat <- data.frame(cbind(Tbs, TempCorr0, TempCorr1, TempCorr2, TempCorr3, damage))
TempCorrDatL <- TempCorrDat %>% pivot_longer(cols = c(TempCorr0, TempCorr1, TempCorr2, TempCorr3),
names_to = c("repair"), values_to = "repair_value")
TempCorrDatL <- subset(TempCorrDatL, Tbs > 6)
TempCorrDatLrepair <- subset(TempCorrDatL, repair!="TempCorr0")
dldat <- data.frame(cbind(data.frame(dl0),data.frame(dl1),data.frame(dl2),data.frame(dl3)))
dldatL <- dldat %>% pivot_longer(cols = c(survival, survival.1, survival.2, survival.3), names_to = c("repair"),
values_to = "repair_value")
dldatKL <- dldat %>% pivot_longer(cols = c(kdot, kdot.1, kdot.2, kdot.3), names_to = c("kdot"),
values_to = "kdot_value")
dldat_comb <- data.frame(cbind(dldatL$repair, dldatL$repair_value, dldatKL$kdot, dldatKL$kdot_value))
colnames(dldat_comb) <- c("repair", "repair_value", "kdot", "kdot_value")
dldat_comb$repair_value <- as.numeric(dldat_comb$repair_value)
dldat_comb$kdot_value <- as.numeric(dldat_comb$kdot_value)
func_seq <- seq(0,100,1)
repair0_func <- kdot0*(kdot0^(3/(func_seq)))
repair1_func <- kdot1*(kdot1^(3/(func_seq)))
repair2_func <- kdot2*(kdot2^(3/(func_seq)))
repair3_func <- kdot3*(kdot3^(3/(func_seq)))
rate_decay <- data.frame(func_seq, repair0_func, repair1_func, repair2_func, repair3_func)
rate_decayL <- rate_decay %>% pivot_longer(cols = c(repair0_func, repair1_func, repair2_func, repair3_func),
names_to = "repair", values_to = "repair_value")
# Plot sensitivity temp~time and display CTmax and z
base_tls_plot <-
ggplot(data = sim, aes(y = ht, x = time2)) +
geom_smooth(method = "lm", formula = y~x, colour = "black", fullrange = F) +
ylim(35,45) + xlim(0,2.6) +
labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
annotate("text", x = 0.3, y = 38, label = bquote(italic(CT[max1h])=="43.8°C"), hjust = 0) +
annotate("text", x = 0.3, y = 36, label = bquote(italic(z)=="-2.4"), hjust = 0) +
theme_classic()
# Damage rate
damage_plot <-
ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = damage)) +
geom_path() +
scale_x_continuous(limits = c(5,45)) +
scale_y_continuous(limits = c(0,0.06)) +
labs(y = bquote(Damage~rate~("%"~min^-1)), x = bquote(Temperature~(degree*C))) +
theme_classic()
# Repair rate
repair_plot <-
ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = repair_value, colour = repair)) +
geom_path(aes(group = rev(repair))) +
scale_x_continuous(limits = c(5,45)) +
scale_y_continuous(limits = c(0,0.06)) +
scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]),
labels = c("None", "Low", "Moderate", "High")) +
labs(y = bquote(Repair~rate~("%"~min^-1)), x = bquote(Temperature~(degree*C))) +
theme_classic()
# Ratio of damage:repair
ratio_plot <-
ggplot(data.frame(TempCorrDatLrepair), aes(x = Tbs, y = (damage/repair_value), colour = repair)) +
geom_path(aes(group = rev(repair))) +
geom_hline(yintercept = 1, linetype = 2) +
scale_x_continuous(limits = c(5,45)) +
scale_y_log10(limits = c(0.005,10000), breaks = c(0.01, 0.1, 1, 10, 100, 1000),
labels = c("0.01", "0.1", "1", "10", "100", "1000")) +
scale_colour_manual(name = "Repair rate", values = c(cbp2[2:4]),
labels = c("Low", "Moderate", "High")) +
annotate("text", x = 32, y = 0.05, label = "Repair \noutweighs \ndamage", hjust = 0) +
annotate("text", x = 12, y = 1000, label = "Damage \noutweighs \nrepair", hjust = 0) +
labs(y = bquote(Damage:repair~ratio), x = bquote(Temperature~(degree*C))) +
theme_classic() + theme(legend.position = "none")
# Net damage rate
net_plot <-
ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = damage-repair_value, colour = repair)) +
geom_path(aes(group = rev(repair))) +
geom_hline(yintercept = 0, linetype = 2) +
scale_x_continuous(limits = c(5,45)) +
scale_y_continuous(limits = c(-0.06,0.1)) +
scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]),
labels = c("None", "Low", "Moderate", "High")) +
annotate("text", x = 32, y = -0.03, label = "Repair \noutweighs \ndamage", hjust = 0) +
annotate("text", x = 12, y = 0.08, label = "Damage \noutweighs \nrepair", hjust = 0) +
labs(y = bquote(Net~damage~rate~("%"~min^-1)), x = bquote(Temperature~(degree*C))) +
theme_classic()
# Plot decline in repair over time due to damage accumulation
decline_plot <-
ggplot(dldatKL, aes(y = kdot_value, x = time, colour = kdot)) +
geom_vline(xintercept = 6350, colour = "red", lty = 2) +
geom_vline(xintercept = 30870, colour = "red", lty = 2) +
geom_vline(xintercept = 36550, colour = "red", lty = 2) +
geom_path(aes(group = kdot)) +
#scale_y_continuous(limits = c(0,100)) +
scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]),
labels = c("None", "Low", "Moderate", "High")) +
labs(x = "Time (days)", y = bquote(Repair~rate~coefficient~(italic(dot(k))))) +
theme_classic() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
# Plot four weeks of body temperature
tb_plot <-
ggplot(data.frame(dl0), aes(y = ta, x = time)) +
geom_hline(yintercept = 28.5, colour = "black", lty = 2) +
geom_vline(xintercept = 6350, colour = "red", lty = 2) +
geom_vline(xintercept = 30870, colour = "red", lty = 2) +
geom_vline(xintercept = 36550, colour = "red", lty = 2) +
geom_line() +
labs(x = "Time (days)", y = bquote(Temperature~(degree*C))) +
scale_y_continuous(limits = c(3,40)) +
annotate("text", x = 12000, y = 38, label = bquote(italic(T)[c]=="28.5°C"), hjust = 0) +
theme_classic() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
# Plot physiological function over the four weeks incorporating repair
func_plot <-
ggplot(dldatL, aes(y = repair_value, x = time, colour = repair)) +
geom_vline(xintercept = 6350, colour = "red", lty = 2) +
geom_vline(xintercept = 30870, colour = "red", lty = 2) +
geom_vline(xintercept = 36550, colour = "red", lty = 2) +
geom_path(aes(group = rev(repair))) +
scale_y_continuous(limits = c(0,100)) +
scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]),
labels = c("None", "Low", "Moderate", "High")) +
labs(x = "Time (days)", y = "Physiological function (%)") +
theme_classic() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
# Plot decline as a function of damage accumulation
decline_func_plot <-
ggplot(rate_decayL, aes(y = repair_value, x = func_seq, colour = repair)) +
geom_path(aes(group = rev(repair))) +
#scale_y_continuous(limits = c(0,100)) +
scale_x_reverse() +
scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]),
labels = c("None", "Low", "Moderate", "High")) +
labs(x = "Physiological function (%)", y = bquote(Repair~rate~coefficient~(italic(dot(k))))) +
theme_classic()
ggarrange(base_tls_plot, damage_plot, repair_plot,
ratio_plot, net_plot, decline_func_plot,
tb_plot, func_plot, decline_plot,
labels = c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)", "(i)"),
align = "hv", ncol = 3, nrow = 3,
common.legend = TRUE, legend = "bottom")
#ggsave("Fig2_R1.pdf", width = 10, height = 10)# Load in climate data and crown dieback data
env <- read.csv("Summer_Climate_Penrith.csv", header = TRUE)
env$date <- dmy(as.character(env$Date))
env.noNA <- env[!is.na(env$Rainfall_mm),] # remove missing values from one column
scaleFactor <- max(env.noNA$Max_T_C) / max(env.noNA$Rainfall_mm)
die <- read.csv("Crown_Dieback.csv", header = TRUE)
die$Date <- dmy(as.character(die$Date))
fig3a <- ggplot(env, aes(x = date)) +
geom_rect(aes(NULL, NULL, xmin = as.Date("2019-12-5"), xmax = as.Date("2020-1-15")), ymin = 45, ymax = 55,
fill = alpha("burlywood", 0.1)) +
geom_rect(aes(NULL, NULL, xmin = as.Date("2020-1-16"), xmax = as.Date("2020-2-10")), ymin = 45, ymax = 55,
fill = alpha("lightblue", 0.1)) +
geom_rect(aes(NULL, NULL, xmin = as.Date("2020-2-11"), xmax = as.Date("2020-4-1")), ymin = 45, ymax = 55,
fill = alpha("seagreen3", 0.05)) +
geom_line(aes(y = Max_T_C), colour = "black", linewidth = 0.5) +
geom_bar(aes(y = Rainfall_mm * scaleFactor), fill = alpha("blue", 0.4),
stat = "identity", width = 0.8, size = 0.25) +
coord_cartesian(xlim=c(as.Date("2019-12-15"), as.Date("2020-03-26")), ylim = c(0, 56)) +
scale_x_date(date_labels = "%d-%b-%y") +
scale_y_continuous(expand = c(0, 0),
limits = c(0, 56),
sec.axis = sec_axis(~./scaleFactor, name="Rainfall (mm)")) +
labs(x = "Date", y = expression("Daily Maximum"~italic("T")[air]~~(degree*C))) +
theme(axis.title.x = element_blank(), axis.text.x = element_blank()) +
annotate("text", x = as.Date("2019-12-29"), y = 51,
label = "Too dry to repair \n heat stress \n damage",
size = 3, colour = "orangered") +
annotate("text", x = as.Date("2020-1-28"), y = 51,
label = "Unfavourable \n conditions \n for repair",
size = 3, colour = "blue3") +
annotate("text", x = as.Date("2020-3-5"), y = 51,
label = "Conditions allow repair \n of damage and canopy \n dieback regeneration",
size = 3, colour = "darkgreen") +
theme_classic()
fig3a
#ggsave("Fig3a_R1.pdf", height = 5, width = 5)
fig3b <-
ggplot(die, aes(x = Date, y = Survival_percent, colour = Species_Code)) +
geom_point(size = 2) +
geom_line(aes(group = Species_Code)) +
coord_cartesian(xlim = c(as.Date("2019-12-15"), as.Date("2020-03-26")), ylim = c(0, 100)) +
scale_x_date(date_labels = "%d-%b-%y") +
scale_colour_manual(values = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3")) +
labs(x="Date", y="Tree crown survival (%)") +
annotate("text", x = as.Date("2020-1-12"), y = 95, label = expression(italic("Acer rubrum")),
parse = T, size = 3, colour = "#E41A1C") + #add text to figure
annotate("text", x = as.Date("2020-3-1"), y = 10, label = expression(italic("Banksia integrifolia")),
parse = T, size = 3, colour = "#377EB8") + #add text to figure
annotate("text", x = as.Date("2020-3-1"), y = 16, label = expression(italic("Liriodendron tulipifera")),
parse = T, size = 3, colour = "#4DAF4A") + #add text to figure
annotate("text", x = as.Date("2020-1-28"), y = 81, label = expression(italic("Syzygium floribundum")),
parse = T, size = 3, colour = "#984EA3") + #add text to figure
theme_classic() +
theme(legend.position = "none")
fig3b
#ggsave("Fig3b_R1.pdf", height = 5, width = 5)# Simulate basic survival curves by life stage
time3 <- seq(1,19,1)
seed <- c(rep(1,8),c(1,0.9,0.7,0.6,0.5,0.45,0.4,0.4,0.4,0.4,0.35))
seedling <- c(rep(1,6),c(0.9,0.8,0.8,0.7,0.6,0.4,0.25,0.2,0.2,0.15,0.1,0.1,0.1))
immature <- c(rep(1,8),c(1,1,0.9,0.9,0.9,0.85,0.85,0.8,0.7,0.65,0.65))
adult <- c(rep(1,8),c(1,1,1,1,1,1,1,1,0.8,0.7,0.7))
demodat_c <- data.frame(cbind(time3, seed, seedling, immature, adult))
demodat2 <- pivot_longer(demodat_c, cols = c(seed, seedling, immature, adult), names_to = "stage")
demodat2$treatment <- "control"
seed <- c(rep(1,7),c(0.9,0.8,0.8,0.6,0.5,0.4,0.3,0.3,0.3,0.3,0.3,0.25)*0.5)
seedling <- c(rep(1,6),c(0.8,0.7,0.7,0.5,0.2,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1)*0.3)
immature <- c(rep(1,8),c(1,0.7,0.6,0.6,0.6,0.55,0.55,0.5,0.45,0.45,0.45)*0.7)
adult <- c(rep(1,10),c(1,0.8,0.8,0.8,0.75,0.7,0.6,0.55,0.55)*0.8)
demodat_h <- data.frame(cbind(time3, seed, seedling, immature, adult))
demodat3 <- pivot_longer(demodat_h, cols = c(seed, seedling, immature, adult), names_to = "stage")
demodat3$treatment <- "heat"
demodat4 <- rbind(demodat2,demodat3)
demodat4$stage <- factor(demodat4$stage, levels = c("seed","seedling","immature","adult"))
fig3d <-
ggplot(demodat4, aes(x = time3, y = value, colour = stage, lty = treatment)) +
geom_line() +
scale_colour_manual(values = cbp2[c(2,3,4,6)], name = "") +
scale_linetype_manual(values = c(1,2), labels = c("Tolerant", "Sensitive"), name = "") +
labs(y = "Survival probability to reach \n next life stage or to reproduce",
x = "Cumulative thermal load") +
theme_classic() + theme(axis.text = element_blank(), legend.position = "bottom") +
guides(colour = guide_legend(nrow = 2, byrow = TRUE), linetype = guide_legend(nrow = 2, byrow = TRUE))
fig3d
#ggsave("Fig3d_R1.pdf", width = 5, height = 5)# Simulate temperature and multi-stressor data
simTemp <- c(30:70)
simTime <- seq(0,2.5,0.061)
simH <- seq(55,35,-0.5) + rnorm(seq(55,35,-0.5),0,0.5)
simA <- seq(50,30,-0.5) + rnorm(seq(50,30,-0.5),0,0.5)
simB <- seq(55,25,-0.74) + rnorm(seq(55,25,-0.74),0,0.5)
simP <- seq(50,20,-0.74) + rnorm(seq(50,20,-0.74),0,0.5)
simAB <- seq(60,28,-0.8) + rnorm(seq(60,28,-0.8),0,0.7)
simTA <- seq(55,29,-0.65) + rnorm(seq(55,29,-0.65),0,0.5)
simTB <- seq(55,25,-0.75) + rnorm(seq(55,25,-0.75),0,0.5)
simTC <- seq(52,12,-1) + rnorm(seq(52,12,-1),0,1)
simdata <- data.frame(simTime, simTemp, simH, simA, simB, simP, simAB, simTA, simTB, simTC)
simdata <- simdata %>% filter(row_number() %% 2 != 1)
simdataA <- simdata %>% pivot_longer(cols = c(simH, simA, simB, simAB), names_to = "sim_type")
simdataA1 <- simdata %>% pivot_longer(cols = c(simH, simP), names_to = "sim_type")
simdataB <- simdata %>% pivot_longer(cols = c(simH, simTA, simTB, simTC), names_to = "sim_type")
sim_mod <- lm(value ~ simTime * sim_type, simdataA1)
simH_mod <- lm(value ~ simTime, subset(simdataA1, sim_type=="simH"))
simP_mod <- lm(value ~ simTime, subset(simdataA1, sim_type=="simP"))
summary(sim_mod)##
## Call:
## lm(formula = value ~ simTime * sim_type, data = simdataA1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9067 -0.2909 -0.1095 0.2975 1.0898
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.1223 0.2187 252.06 <2e-16 ***
## simTime -8.1656 0.1553 -52.59 <2e-16 ***
## sim_typesimP -5.2454 0.3093 -16.96 <2e-16 ***
## simTime:sim_typesimP -3.8012 0.2196 -17.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4885 on 36 degrees of freedom
## Multiple R-squared: 0.9972, Adjusted R-squared: 0.997
## F-statistic: 4266 on 3 and 36 DF, p-value: < 2.2e-16summary(simH_mod)##
## Call:
## lm(formula = value ~ simTime, data = subset(simdataA1, sim_type ==
## "simH"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.89227 -0.23981 -0.02263 0.28588 0.84318
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.1223 0.1995 276.34 <2e-16 ***
## simTime -8.1656 0.1416 -57.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4456 on 18 degrees of freedom
## Multiple R-squared: 0.9946, Adjusted R-squared: 0.9943
## F-statistic: 3323 on 1 and 18 DF, p-value: < 2.2e-16summary(simP_mod)##
## Call:
## lm(formula = value ~ simTime, data = subset(simdataA1, sim_type ==
## "simP"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9067 -0.2947 -0.1378 0.3041 1.0898
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.8769 0.2363 211.04 <2e-16 ***
## simTime -11.9668 0.1678 -71.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.528 on 18 degrees of freedom
## Multiple R-squared: 0.9965, Adjusted R-squared: 0.9963
## F-statistic: 5085 on 1 and 18 DF, p-value: < 2.2e-16fig3e <-
ggplot(data = simdataA, aes(x = simTime, y = value, colour = sim_type)) +
geom_smooth(aes(y = value), method = "lm", formula = y~x, fullrange = T, se = F) +
labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
coord_cartesian(ylim = c(41,60)) +
scale_x_continuous(expand = c(0, 0), limits = c(0,2.5)) +
scale_colour_manual(values = c(cbp2[8],cbp2[3],cbp2[4],cbp2[1]),
labels = c(bquote(Stressor~altering~italic(CT[max1h])),
bquote(Stressor~altering~italic(CT[max1h])~and~italic(z)),
bquote(Stressor~altering~italic(z)),
"Heat stress only"), name = "") +
theme_classic() + theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 2, byrow = TRUE), linetype = guide_legend(nrow = 2, byrow = TRUE))
fig3e
#ggsave("Fig3e_R1.pdf", height = 5, width = 5)
fig3f <-
ggplot(data = simdataB, aes(x = simTime, y = value, colour = sim_type)) +
geom_smooth(aes(y = value, linetype = sim_type), method = "lm", formula = y~x, fullrange = T, se = F) +
labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
coord_cartesian(ylim = c(28,60)) + scale_x_continuous(expand = c(0, 0), limits = c(0,2.5)) +
scale_colour_manual(values = c(cbp2[1],cbp2[2],cbp2[7],cbp2[7]),
labels = c("Heat stress only",
"Heat stress + Stressor A",
"Heat stress + Stressor B",
"Heat stress + Stressors A + B"), name = "") +
scale_linetype_manual(values = c(1,3,3,1),
labels = c("Heat stress only",
"Heat stress + Stressor A",
"Heat stress + Stressor B",
"Heat stress + Stressors A + B"), name = "") +
theme_classic() + theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 2, byrow = TRUE), linetype = guide_legend(nrow = 2, byrow = TRUE))
fig3f
#ggsave("Fig3f_R1.pdf", height = 5, width = 5)# Generate tolerance landscape for parameterising models
# Drosophila suzukii
data <- read.csv("Dsuzukii_data_Orsted2024.csv")
head(data)## id temp time sex rep prod dead t_coma
## 1 38-300-F-1 38 48.30624 FALSE 1 0 1 14.50
## 2 38-300-F-2 38 48.30624 FALSE 2 0 1 20.80
## 3 38-300-F-3 38 48.30624 FALSE 3 0 1 13.42
## 4 38-300-F-4 38 48.30624 FALSE 4 0 1 16.25
## 5 38-300-F-5 38 48.30624 FALSE 5 0 1 12.50
## 6 38-300-F-6 38 48.30624 FALSE 6 0 1 7.50tl <- tolerance.landscape(data$temp, data$time)
Tb <- data$temp
time <- data$time
surv.time <- tl$S[ ,2]
surv.pct <- seq(100, 0, length.out = length(c(0, surv.time)))
par(mfrow = c(1, 1), mar = c(4.5, 4, 1, 1))
plot(c(0, surv.time), surv.pct, type = 'l', ylab = 'survival, %', xlab = 'time')
Tb.mn <- tl$ta.mn
z <- 3.27 # tl$z
slope <- -1/z
intercept <- 11.092 # -tl$ctmax * slope for hours. # Set to 11.092 to replicate ctmax.4h value from Ørsted et al 2024
plot(seq(30, 41), 10 ^ (seq(30, 41) * slope + intercept),
ylab = 'survial time, hours', xlab = 'body temperature, C', pch = 16, type = 'b')
ctmax.1h <- (log10(1) - intercept) / slope # temperature causing 50% mortality in 1 hour
ctmax.4h <- (log10(4) - intercept) / slope
ctmax.24h <- (log10(24) - intercept) / slope # temperature causing 50% mortality in 1 day
R0 <- 1
k <- log(10) / z
# Set Arrhenius model and repair rate parameters
TA <- 14065 # Arrhenius temperature, K
TL <- 10.5 + 273.15 # Arrhenius temperature lower threshold, K
TH <- 28.5 + 273.15 # Arrhenius temperature upper threshold, K
TAL <- 50000 # Arrhenius temperature lower, K
TAH <- 100000 # Arrhenius temperature upper, K
TREF <- 20 + 273.15 # reference temperature, K
Tbs <- seq(0, 50) # sequence of temperatures over which to model, deg C
kdot_ref <- 0.00 # repair rate % per minute at ref temp 20 deg C
# Plot the thermal response curve for repair
repair1 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
par(mfrow = c(1, 1), mar = c(4.5, 4, 1, 1))
plot(Tbs, repair1, type = 'l', ylab = 'repair, % per min', xlab = 'body temperature')
abline(v = 28.5, col = 2)
#Tb <- seq(30, 37,1/(60 * 24))
#Tb <- rep(37, 60 * 24)
Tb <- rep(ctmax.1h, 60 * 24)
shift <- 10 ^ ((Tb.mn - Tb) / z)
# minutes to shift survival curve by at each minute given temperature at that minute
time.rel <- 0
result <- matrix(data = 0, nrow = length(Tb) - 1, ncol = 2)
repair <- ArrFunc5(Tb, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
for(i in 1:(length(Tb) - 1)){
#for(i in 130:150 - 1){
# adjust time for the survival probability vs time relationship to account
# for current temperature and then get % alive after the current time increment
alive <- try(approx(c(0, shift[i] * surv.time), surv.pct,
xout = time.rel + 1, ties = "ordered")$y, silent = TRUE)
alive <- min(100, alive + repair[i]) # add repair, but still catch NA
if(is.na(alive)){ # cap at 100% surviving
alive <- 100
break
}
result[i, 1]<- alive # save current survival fraction
# now interpolate to get the time that goes with the survival %
time.rel <- try(approx(surv.pct, c(0, shift[i + 1] * surv.time),
xout = alive)$y, silent = TRUE)
result[i, 2] <- time.rel # save time since start
}
plot(result[, 2], result[, 1], type = 'l', ylab = "% surviving", xlab = "time, min")
# Integrate damage (experimental end point of death/coma gives lethal dose)
L <- matrix(0, nrow = length(Tb))
for(i in 2:(length(Tb))){
if(Tb[i] > ctmax.24h){
dLdt <- exp(k * (Tb[i] - ctmax.24h) - 1)
}else{
dLdt <- 0
}
L[i] <- L[i-1] + dLdt
}
Ld <- L[which(result[, 1]==0)[1]]
Ld## [1] 1748.163Spatial distribution simulation
# day of year to start simulation
DOY <- 1
sim.length <- 7 # n days to simulate
# read in air temps and filter for month, ctmax
tair.120cm <- brick("TA120cm_2017.nc")## Error in charToDate(x) :
## character string is not in a standard unambiguous formattair.120cm.sub <- tair.120cm[[(DOY * 24 + 1):(DOY * 24 + sim.length * 24)]]
tair.ctmax <- tair.120cm.sub / 10
tair.ctmax[tair.ctmax < ctmax.1h] <- 0
tair.ctmax[tair.ctmax > ctmax.1h] <- 1
sum.sub <- calc(tair.ctmax, function(x){sum(x)})
max.sub <- calc(tair.120cm.sub, function(x){max(x)})
ct.range <- sum.sub
ct.range[ct.range == 0] <- -1
ct.range[ct.range >= 0] <- 0
ct.range[ct.range == -1] <- 1
# create grid of points
nc <- nc_open("TA120cm_2017_time.nc")
lon <- ncvar_get(nc, "longitude") # extract longitude values
lat <- ncvar_get(nc, "latitude") # extract latitude values
sites <- expand.grid(lon, lat) # 3618 sites to iterate over each day
# Loop across sites and calculate damage (first with no repair)
# for dynamic tolerance landscapes (Rezende approach) and dynamic ctmax (Jørgensen approach) using the
# Drosophila suzukii female fecundity parameterised data plotted over Australia.
# Note: Simulation may take well over 1 hour to run for each level of repair (0 or 0.003)
kdot_ref <- 0 ### change to 0.02 to simulate repair, or 0 for no repair ###
all.results <- matrix(nrow = nrow(sites), ncol = 4)
par(mfrow = c(1, 1), mar = c(4.5, 4, 1, 1))
plot(ct.range, zlim = c(0.5, 1), col = 'grey', main = "range based on 120 cm air temperature, deg C",
ylab = "latitude", xlab = "longitude")
ozmap(x = "country", add = TRUE)
Code for running spatial distribution simulation (not run in markdown)
# parameters referring to survival or surv can be any sublethal measure, this is a generic function term
kdot_ref <- 0
for(ii in 1:nrow(sites)){ # start site loop
loc <- c(sites[ii, 1], sites[ii, 2])
tair.120cm_loc <- extract.nc("TA120cm_2017_time.nc", loc)[(DOY * 24 + 1):(DOY * 24 + sim.length * 24), ]
if(!is.na(tair.120cm_loc$data[1])){ # check if land or sea
if(is.na(all.results[ii, 1])){
variable.temp <- tair.120cm_loc$data / 10
if(max(variable.temp) < tl$ctmax){ # don't do sites that won't survive one minute
time.min <- spline(variable.temp, n = length(variable.temp - 1) * 60)$x - 1
Tb.min <- spline(variable.temp, n = length(variable.temp - 1) * 60)$y
# Daily survival probability
# dynamic.landscape function (function returns NA for survival probability < 0)
# Loop every 1440 min (= 24 h * 60 min)
day <- rep(1:sim.length, each = 1440)
surv.out <- seq(100, 0, length.out = length(c(0, surv.time))) # vector to hold survival
# plot(surv.time, surv.pct[-1], ylab = 'surviving %', xlab = 'time, min', type = 'l')
final <- matrix(data = 100, nrow = sim.length, ncol = 1)
final2 <- final
alive <- 100
last.alive <- alive
damage <- 0
result2 <- matrix(data = 0, nrow = length(tair.120cm_loc$data) * 60, ncol = 2)
for(xx in 1:sim.length){ # loop through days
Tb <- Tb.min[xx == day]
shift <- 10 ^ ((Tb.mn - Tb) / z) # minutes to shift survival curve by at each minute
time.rel <- 0
result <- matrix(data = 0, nrow = length(Tb) - 1, ncol = 2)
repair <- ArrFunc5(Tb, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
#plot(seq(0, 1440-1), Tb, ylab = 'body temperature', xlab = 'time, min', type = 'l')
#plot(seq(0, 1440-1), repair, ylab = '% repair / min', xlab = 'time, min', type = 'l')
#plot(seq(0, 1440-1), shift, ylab = 'survival shift, min', xlab = 'time, min', type = 'l')
# count <- 0
if(xx == 1){
starti = 2
result[1, 1] <- 100
}else{
starti = 1
}
for(i in starti:(length(Tb) - 1)){
alive <- try(approx(c(0, shift[i] * surv.time),
surv.pct, xout = time.rel + 1, ties = "ordered")$y,
silent = TRUE)
if(is.na(alive)){
if(last.alive != 100){
alive <- 0
}else{
alive <- 100
}
}
if(alive == 0){
break
}
alive <- min(100, alive + repair[i]) # add repair, but still catch NA
# modify kdot_ref based on alive. If alive > 99, kdot_ref stays the same.
# If alive < 99, kdot_ref reduces as a function of alive
kdot_ref1 <- ifelse(alive < 99, kdot_ref*(kdot_ref^(3/(alive))), kdot_ref)
alive <- alive + ArrFunc5(Tb[i], TA, TAL, TAH, TL, TH, TREF, kdot_ref1)
alive <- ifelse(alive > 100, 100, alive)
last.alive <- alive
result2[(xx - 1) * 24 * 60 + i, ] <- c(i, alive)
result[i, 1] <- alive
time.rel <- try(approx(surv.pct, c(0, shift[i + 1] * surv.time), xout = alive)$y, silent = TRUE)
result[i, 2] <- time.rel
}
cat(xx, '\n')
final[xx] <- min(result[, 1], na.rm = TRUE)
# Jørgensen et al. approach - adapted from equation 2: https://doi.org/10.1038/s41598-021-92004-6
L <- matrix(Ld, nrow = length(Tb))
L[1] <- damage
for(i in 2:(length(Tb))){
if(Tb[i] > ctmax.24h){
dLdt <- exp(k * (Tb[i] - ctmax.24h) - 1)
}else{
dLdt <- 0
}
L[i] <- max(0, L[i-1] + dLdt - repair[i] / 100 * Ld)
if(L[i] > Ld){ # dead
L[i] <- Ld
break
}
}
#points(L/Ld * 100, col = 'purple', type = 'l')
# work out damage end point for carry over to next day
damage.end <- tail(L, 1)
if(damage.end > Ld){
damage <- Ld
}
final2[xx] <- max(L, na.rm = TRUE)
} # end of day loop
#plot(100 - result2[, 2], type = 'l')
#plot(100 - final, type = 'h')
#points(final2 / Ld * 100, type = 'p', pch = 16, col = 2, ylim = c(0, 100))
#plot(100-final, final2 / Ld * 100, pch = 16, ylim = c(0, 100), xlim = c(0, 100), ylab = 'Jorgensen', xlab = 'Rezende')
#abline(0, 1)
points(loc[1], loc[2], cex = sum(final[,1])/30/100, pch = 16, col = 'blue')
if(min(final[,1]) / 75 == 0){
points(loc[1], loc[2], pch = 4)
}else{
points(loc[1], loc[2], cex = min(final[,1]) / 75, pch = 16, col = 'blue')
}
# write.table(c(loc, min(final[,1])), file = "Rezende.out.rev0.csv", sep = ",",
# col.names = F, qmethod = "double", append = T)
# write.table(c(loc, max(final2[,1])/Ld * 100), file = "Jorgensen.out.rev0.csv", sep = ",",
# col.names = F, qmethod = "double", append = T)
all.results[ii, ] <- c(loc, min(final[, 1]), max(final2[, 1]) / Ld * 100)
cat('site', ii, '\n')
}else{
points(loc[1], loc[2], pch = 4)
}
}
}
}
all.results2 <- as.data.frame(na.omit(all.results))
colnames(all.results2) <- c('lon', 'lat', 'surv', 'dose')
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE)
points(all.results2[, 1:2], cex = 1 - all.results2$dose / 75, pch = 16, col = 'blue')
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE)
points(all.results2[, 1:2], cex = all.results2$surv / 120, pch = 16, col = 'blue')
all.results.no.repair <- all.results2
#write.csv(all.results.no.repair, 'dsuz.results.no.repair0.csv')
#write.csv(all.results.with.repair, 'dsuz.results.repair0.02.csv')dsuz.results.no.repair0 <- read.csv("dsuz.results.no.repair.pro0.csv")[, -1]
dsuz.results.repair2 <- read.csv("dsuz.results.repair.pro0.02.csv")[, -1]
rezende.no.repair <- read.csv("Rezende.out.pro0.csv", header = F, col.names = c("Var","Value"))
jorgensen.no.repair <- read.csv("Jorgensen.out.pro0.csv", header = F, col.names = c("Var","Value"))
head(rezende.no.repair)## Var Value
## 1 1 132.11900
## 2 2 -11.40100
## 3 3 99.41495
## 4 1 142.31900
## 5 2 -11.40100
## 6 3 99.43911head(jorgensen.no.repair)## Var Value
## 1 1 132.119
## 2 2 -11.401
## 3 3 0.000
## 4 1 142.319
## 5 2 -11.401
## 6 3 0.000rezende.no.repair_wide <- pivot_wider(rezende.no.repair, names_from = Var, values_from = Value)
rezende.no.repair_wide <- cbind(data.frame(rezende.no.repair_wide$`1`, rezende.no.repair_wide$`2`, rezende.no.repair_wide$`3`))
names(rezende.no.repair_wide) <- c("lat", "lon", "r_surv")
head(rezende.no.repair_wide)## lat lon r_surv
## 1 132.119 -11.401 99.41495
## 2 142.319 -11.401 99.43911
## 3 132.719 -12.001 77.52720
## 4 133.319 -12.001 49.95959
## 5 133.919 -12.001 66.39607
## 6 142.319 -12.001 88.41603jorgensen.no.repair_wide <- pivot_wider(jorgensen.no.repair, names_from = Var, values_from = Value)
jorgensen.no.repair_wide <- cbind(data.frame(jorgensen.no.repair_wide$`1`, jorgensen.no.repair_wide$`2`, jorgensen.no.repair_wide$`3`))
names(jorgensen.no.repair_wide) <- c("lat", "lon", "j_surv")
head(jorgensen.no.repair_wide)## lat lon j_surv
## 1 132.119 -11.401 0.000000
## 2 142.319 -11.401 0.000000
## 3 132.719 -12.001 16.770258
## 4 133.319 -12.001 42.385432
## 5 133.919 -12.001 27.874079
## 6 142.319 -12.001 8.570505no.repair.dat <- as.data.frame(cbind(rezende.no.repair_wide, (100-jorgensen.no.repair_wide[,3])))
names(no.repair.dat) <- c("lat", "lon", "r_surv", "j_surv")
no.repair.dat$diff <- no.repair.dat$r_surv - no.repair.dat$j_surv
no.repair.dat$diff.pos <- ifelse(no.repair.dat$diff < 0, NA, no.repair.dat$diff)
no.repair.dat$diff.neg <- ifelse(no.repair.dat$diff > 0, NA, no.repair.dat$diff)
tair.brks <- seq(20, 50, 5) # breaks for the colour scheme
tair.col <- colorRampPalette(c("blue", "yellow", "red"))(length(tair.brks) - 1) # colour# Multi-panel plot
par(mfrow = c(2,2), mar = c(3,1,3,1), oma = c(1,1,3,1))
# Maximum air temp over the simulation
plot(max.sub / 10, breaks = tair.brks, col = tair.col, ylab = 'latitude', xlab = 'longitude', axes = F, box = F,
main = bquote(atop(bold(a)~"Maximum air temperature"~(120~cm~","~degree*C),
Productivity~probability~("%")~italic(CT["max1h"])~model)), alpha = 0.3)
plot(max.sub / 10, zlim = c(20, ctmax.1h), breaks = tair.brks, col = tair.col,
ylab = 'latitude', xlab = 'longitude', add = T, legend = F)
ozmap(x = "country", add = TRUE)
# Rezende model
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
main = bquote(atop(bold(b)~Productivity~probability~("%"),~without~repair~model)))
points(dsuz.results.no.repair0[, 1:2], cex = dsuz.results.no.repair0$surv/120, pch = 16, col = cbp2[6])
ozmap(x = "country", add = TRUE)
legend(x = 124, y = -32, legend = paste(seq(20, 100, 20), "%"), pt.cex = seq(20, 100, 20)/75, pch = 16, col = cbp2[6], bty = 'n')
# Jørgensen model
# plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
# main = bquote(atop(bold(b)~Survival~probability~("%")~without~repair,Jørgensen~italic(et~al.)~dynamic~italic(CT[max])~model)))
# points(no.repair.dat[, 1:2], cex = no.repair.dat$j_surv/120, pch = 16, col = cbp2[3])
# ozmap(x = "country", add = TRUE)
# legend(x = 126, y = -32, legend = paste(seq(20, 100, 20), "%"), pt.cex = seq(20, 100, 20)/75, pch = 16, col = cbp2[3], bty = 'n')
# Rezende model with repair at kdot=0.02
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
main = bquote(atop(bold(c)~Productivity~probability~("%"),~with~high~repair~rate~model)))
points(dsuz.results.repair2[, 1:2], cex = dsuz.results.repair2$surv / 120, pch = 16, col = cbp2[4])
ozmap(x = "country", add = TRUE)
legend(x = 124, y = -32, legend = paste(seq(20, 100, 20), "%"), pt.cex = seq(20, 100, 20)/75, pch = 16, col = cbp2[4], bty = 'n')
# Diff between no repair and repair at 0.02
dsuz.results.repair2$diff <- dsuz.results.repair2$surv - dsuz.results.no.repair0$surv
range(dsuz.results.repair2$diff)## [1] -10.59975 10.67527plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
main = bquote(atop(bold(d)~Productivity~probability~difference~("%"),~with~repair~"-"~without~repair~models)))
points(dsuz.results.repair2[, 1:2], cex = dsuz.results.repair2$diff/15, pch = 16, col = cbp2[8])
ozmap(x = "country", add = TRUE)
legend(x = 124, y = -32, legend = paste(seq(5, 25, 5), "%"), pt.cex = seq(5, 25, 3)/15, pch = 16, col = cbp2[8], bty = 'n')
# Save figure as approximately 8 x 8 inch